#loading packages
library('FactoMineR')
library('factoextra')
library("gplots")
library("ggrepel")
library("dplyr")
library("tidyverse")
library("ggplot2")
library("ggthemes")
library("forcats")
library("reshape2")

#datasets

#Chapter 2
load(file="Data/Existent Data on Religion/Treated Data/Final Data/ReligiousAccommodationEurope.rdata") #Dataset used for the MCA - Country MCA

load(file="Data/Existent Data on Religion/Treated Data/Final Data/accommeans.rdata") #Dataset used for the Score analysis - Country MCA

#MCA for Countries - ReligiousAccomodationEurope dataset

##load dataset
load(file="Data/Existent Data on Religion/Treated Data/Final Data/ReligiousAccommodationEurope.rdata") #Dataset used for the MCA - Country MCA

##MCA
mca <- MCA(Accom_active, graph = FALSE)
print(mca)

##Variable Coordinates Contribution
var <- get_mca_var(mca)
dt <-var$coord
dt <- data.frame(dt[1:36, 1:2])

##Adding Level of Restriction and Contribution
Restriction <- c(0, 1, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 0, 2, 3, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 0, 2, 0, 2, 0, 1, 2, 3)
dt$Restrictions <- Restriction
dt$Mean <- round(rowMeans(dt[,c('Dim.1', 'Dim.2')], na.rm=TRUE), 2)
dt$Size <- abs(dt$Mean)
size <- c(0, 0, 0, 3, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0, 0, 0, 1, 0, 0, 1, 1)
dt$size <- size

rownames(dt) <- c("Rites of Passage", "Rites of Passage ", "Diet Laws", "Diet Laws ", "Diet Laws  ", "Burial", "Burial ", "Burial  ", "Rel. Clothing", "Rel. Clothing ", "Rel. Clothing  ", "Rel. Symbols", "Rel. Symbols ", "Public Holidays", "Public Holidays ", "Public Holidays  ", "Private Holidays", "Private Holidays ", "Clergy Jails", "Clergy Jails ", "Clergy Jails  ", "Clergy Military Base", "Clergy Military Base ", "Clergy Military Base  ", "Clergy Military Base   ", "Clergy Hospitals", "Clergy Hospitals ", "Clergy Hospitals  ", "Public Preaching", "Public Preaching ", "Proselytizing", "Proselytizing ", "Religious Ed. Mandatory", "Religious Ed. Mandatory ", "Religious Ed. Mandatory  ", "Religious Ed. Mandatory   ")

##Visualising Variable Categories

ggplot(dt, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color = factor(Restrictions), size=Size)) +
  labs(x = "Privatisation of Religion", y = "Regulation of Religion") +
  geom_hline(yintercept=0, color="grey") +
  geom_vline(xintercept=0, color="grey") +
  geom_text_repel(label=rownames(dt), size=5)+
  theme_clean() +
  labs(title = "Fig.1 Religious Accommodation Variable Categories", caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration") +
  theme(title= element_text(size = 15), text= element_text(size= 8), legend.position = c(0.2, 0.25), legend.justification = c("right", "top"),  legend.text = element_text(size=15), axis.text=element_text(size=12), axis.title=element_text(size=14)) +
  guides(color = guide_legend(override.aes = list(size=5))) +
  labs(color="Degree of Intervention") +
  scale_color_manual(labels = c("No", "Little", "Somewhat", "High"), values=c("#488f31", "#9fc08f", "#fda0c3", "#fa3497"))+
  guides(size=FALSE)

##Extracting Individuals Data for Country Categorisation
df <-mca$ind$coord
df <- data.frame(df[1:216, 1:2])

##Selecting and summarising per countries
df_test<- data.frame(df[1:8, 1:2])
df_test <- summarise_each(df_test, funs(mean))
rownames(df_test) <- c("Austria")

Belgium <- data.frame(df[9:16, 1:2])
Belgium <- summarise_each(Belgium, funs(mean))
rownames(Belgium) <- c("Belgium")
Clusters <- dplyr::union(df_test, Belgium)

Bulgaria <- data.frame(df[17:24, 1:2])
Bulgaria <-summarise_each(Bulgaria, funs(mean))
rownames(Bulgaria) <- c("Bulgaria")
Clusters <- dplyr::union(Clusters, Bulgaria)

Croatia <- data.frame(df[25:32, 1:2])
Croatia <- summarise_each(Croatia, funs(mean))
rownames(Croatia) <- c("Croatia")
Clusters <- dplyr::union(Clusters, Croatia)

Czech <- data.frame(df[33:40, 1:2])
Czech <- summarise_each(Czech, funs(mean))
rownames(Czech) <- c("Czech")
Clusters <- dplyr::union(Clusters, Czech)

Denmark <- data.frame(df[41:48, 1:2])
Denmark <- summarise_each(Denmark, funs(mean))
rownames(Denmark) <- c("Denmark")
Clusters <- dplyr::union(Clusters, Denmark)

Estonia <- data.frame(df[49:56, 1:2])
Estonia <- summarise_each(Estonia, funs(mean))
rownames(Estonia) <- c("Estonia")
Clusters <- dplyr::union(Clusters, Estonia)

Finland <- data.frame(df[57:64, 1:2])
Finland <- summarise_each(Finland, funs(mean))
rownames(Finland) <- c("Finland")
Clusters <- dplyr::union(Clusters, Finland)

France <- data.frame(df[65:72, 1:2])
France <- summarise_each(France, funs(mean))
rownames(France) <- c("France")
Clusters <- dplyr::union(Clusters, France)

Germany <- data.frame(df[73:80, 1:2])
Germany <- summarise_each(Germany, funs(mean))
rownames(Germany) <- c("Germany")
Clusters <- dplyr::union(Clusters, Germany)

Greece <- data.frame(df[81:88, 1:2])
Greece <- summarise_each(Greece, funs(mean))
rownames(Greece) <- c("Greece")
Clusters <- dplyr::union(Clusters, Greece)

Hungary <- data.frame(df[89:96, 1:2])
Hungary <- summarise_each(Hungary, funs(mean))
rownames(Hungary) <- c("Hungary")
Clusters <- dplyr::union(Clusters, Hungary)

Ireland <- data.frame(df[97:104, 1:2])
Ireland <- summarise_each(Ireland, funs(mean))
rownames(Ireland) <- c("Ireland")
Clusters <- dplyr::union(Clusters, Ireland)

Italy <- data.frame(df[105:112, 1:2])
Italy <- summarise_each(Italy, funs(mean))
rownames(Italy) <- c("Italy")
Clusters <- dplyr::union(Clusters, Italy)

Latvia <- data.frame(df[113:120, 1:2])
Latvia <- summarise_each(Latvia, funs(mean))
rownames(Latvia) <- c("Latvia")
Clusters <- dplyr::union(Clusters, Latvia)

Lithuania <- data.frame(df[121:128, 1:2])
Lithuania <- summarise_each(Lithuania, funs(mean))
rownames(Lithuania) <- c("Lithuania")
Clusters <- dplyr::union(Clusters, Lithuania)

Luxembourg <- data.frame(df[129:136, 1:2])
Luxembourg <- summarise_each(Luxembourg, funs(mean))
rownames(Luxembourg) <- c("Luxembourg")
Clusters <- dplyr::union(Clusters, Luxembourg)

Malta <- data.frame(df[137:144, 1:2])
Malta <- summarise_each(Malta, funs(mean))
rownames(Malta) <- c("Malta")
Clusters <- dplyr::union(Clusters, Malta)

Netherlands <- data.frame(df[145:152, 1:2])
Netherlands <- summarise_each(Netherlands, funs(mean))
rownames(Netherlands) <- c("Netherlands")
Clusters <- dplyr::union(Clusters, Netherlands)

Poland <- data.frame(df[153:160, 1:2])
Poland <- summarise_each(Poland, funs(mean))
rownames(Poland) <- c("Poland")
Clusters <- dplyr::union(Clusters, Poland)

Portugal <- data.frame(df[161:168, 1:2])
Portugal <- summarise_each(Portugal, funs(mean))
rownames(Portugal) <- c("Portugal")
Clusters <- dplyr::union(Clusters, Portugal)

Romania <- data.frame(df[169:176, 1:2])
Romania <- summarise_each(Romania, funs(mean))
rownames(Romania) <- c("Romania")
Clusters <- dplyr::union(Clusters, Romania)

Slovakia <- data.frame(df[177:184, 1:2])
Slovakia <- summarise_each(Slovakia, funs(mean))
rownames(Slovakia) <- c("Slovakia")
Clusters <- dplyr::union(Clusters, Slovakia)

Slovenia <- data.frame(df[185:192, 1:2])
Slovenia <- summarise_each(Slovenia, funs(mean))
rownames(Slovenia) <- c("Slovenia")
Clusters <- dplyr::union(Clusters, Slovenia)

Spain <- data.frame(df[193:200, 1:2])
Spain <- summarise_each(Spain, funs(mean))
rownames(Spain) <- c("Spain")
Clusters <- dplyr::union(Clusters, Spain)

Sweden <- data.frame(df[201:208, 1:2])
Sweden <- summarise_each(Sweden, funs(mean))
rownames(Sweden) <- c("Sweden")
Clusters <- dplyr::union(Clusters, Sweden)

UK <- data.frame(df[209:216, 1:2])
UK <- summarise_each(UK, funs(mean))
rownames(UK) <- c("UK")
Clusters <- dplyr::union(Clusters, UK)

##Visualising Country Placements Biplot

ggplot(Clusters, aes(x=Dim.1, y=Dim.2)) +
  geom_point(color="#e43085", size=3) +
  labs(x = "Privatisation of Religion", y = "Regulation of Religion") +
  geom_hline(yintercept=0, color="grey") +
  geom_vline(xintercept=0, color="grey") + 
  geom_text_repel(
    label=rownames(Clusters), size=5
  ) +
  theme_clean() +
  labs(title = "Fig.3 Religious Accommodation", subtitle = "Average country placement according to Privatisation of Religion and Regulation of Religion.", caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration") +
  theme(title= element_text(size = 15), text= element_text(size= 8), axis.text = element_text(size=12), axis.title = element_text(size=14))

##Adding all values for ranking
load("Data/Existent Data on Religion/Treated Data/Final Data/accommeans.rdata")
rel <- round(rel, digits=2)

names(rel) <- c("Rites of Passage","Dietary Laws", "Burial","Religious Clothing", "Religious Symbols","Public Holidays", "Private Holidays", "Clergy Jails","Clergy Military Base", "Clergy Hospitals", "Public Preaching","Proselytizing","Mandatory Religious Education")

rel <- rel[,c("Rites of Passage","Dietary Laws", "Burial","Religious Clothing", "Religious Symbols","Public Holidays", "Private Holidays", "Clergy Jails","Clergy Military Base", "Clergy Hospitals", "Public Preaching","Proselytizing","Mandatory Religious Education")]

rel$Sum <- rowSums(rel)
rel <- arrange(rel, desc(Sum))
rel = subset(rel, select = -c(Sum) )

rel$Country <- rownames(rel)

rel <- rel[,c("Country",names(sort(colSums(rel[,c(1:ncol(rel)-1)]), decreasing = T)))]

ord <- rel$Country
ord <- ord[length(ord):1]

rel1 <- melt(rel)

rel1 <- rel1 %>%
  mutate(Country = fct_relevel(Country, ord))

rel1$value <- round(rel1$value,1)

##Heatmap
rel1 %>%
  ggplot( aes(variable, Country, fill= value))+
  geom_tile(color = "white",
            lwd = 0.5,
            linetype = 1) +
  geom_text(aes(label = value), color = "white", size = 5)+
  xlab("Variable") +
  ylab("Country") +
  ggtitle("Fig.4 Average Country Scores in\nOriginal Datasets")+
  labs(fill = "Restrictions") +
  scale_fill_gradient(low="white",high="#fa3497")+
  theme_clean() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 60, vjust = 0, hjust=0, size=12),
        plot.margin=unit(c(0.5,1,0.5,0),"cm"),
        plot.background = element_rect(color=NA), title= element_text(size = 13), text= element_text(size= 13), axis.text.y = element_text(size=12))+
  xlab("")+
  ylab("")+
  scale_x_discrete(position = "top")+
  labs(caption= "Source: Religion and the State Round 3, Global Restrictions on Religion. Own Elaboration")